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We present a novel potential model for calculating the interaction between a molecule and a single 
graphene sheet. The dispersion/repulsion, induction, dipole-quadrupole, quadrupole-quadrupole 
interactions between a fluid molecule and a graphene sheet are described by integrated functions that 
are only dependent on the separation between the molecule and the graphene along its planar normal. 
The derived potential functions are in excellent agreement with the computationally demanding 
atom-explicit summation method. Typical errors of the integrated potential are less than 2% in 
the energy minimum compared with the exact atom-explicit summation. To examine the practical 
effectiveness of the newly developed functions, Monte Carlo simulations were performed to model 
the adsorption of two representative gases in graphene sheets using both the integrated and atom- 
explicit potentials. The integrated potential results in same adsorption isotherms and density profiles 
for the adsorbed phase while it only requires negligible computing time compared with that using 
the atom-explicit method. The newly developed potential functions provide a simple and accurate 
approach to calculating the physical interaction between molecules and graphene sheets. 



I. INTRODUCTION 

Graphene, or individual layers of graphite in which 
each carbon atom is bonded to three other carbon atoms, 
has been the subject of considerable research interest re- 
cently Being a representative two dimensional crys- 
tal, graphene has many peculiar properties such as low 
dimensionality, surface homogeneity, structure stability, 
conductivity, and charge transport ability. It is deemed 
as one of the most promising materials for future applica- 
tions in various fields such as electronics, novel materials, 
sensors, biodevices etc. [Il-tl3|. 

In the investigations of fundamental properties of 
graphene, it is often of great interest to understand the 
interaction of molecules with graphene on a molecule 
level. For example, graphene sheets have been used as 
high-sensitivity sensors [8| or background membranes 14 1 
in studying the behavior of individual molecules. In such 
applications, interaction between the target molecule and 
the graphene is one of the key properties that needs 
to be understood in order to design effective devices. 
Therefore, a simple and reliable approach to estimating 
the graphene-fluid interaction is essential. Furthermore, 
computer simulation, which is becoming a useful part- 
ner of experiments in studying fluid-solid interactions, 
requires realistic and computationally affordable method 
with sufficient accuracy to calculate the molecule-solid in- 
teractions, such as that between molecules and graphene 
sheets. However, to our knowledge, a simple and accurate 
function describing the interaction between a molecule 
and a single graphene sheet is still not available. In this 
paper, we attempt to develop a novel potential to solve 
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this problem. 

There are two general approaches to calculating the 
interactions between a molecule and a crystal solid like 
graphene. One can use an atom-explicit potential to de- 
scribe the interaction between each atom in the molecule 
and each atom in crystal, calculating the total potential 
energy by summing up all pairs of atoms. Dispersion, 
repulsion, electrostatic and induction interactions can all 
be computed in this way. But such an approach is com- 
putationally costly, especially for large systems involving 
significant number of solid atoms. Alternatively, one can 
develop integrated potential functions that accounts for 
all the solid atoms in an effective way, making use of 
the periodicity and homogeneity in the distribution of 
atoms in the crystal. This has been a popular and useful 
approach for modeling solid-fluid interactions for many 
systems, including graphite and metals [15j . Integrated 
potentials are easy to program and very efficient compu- 
tationally. 

In this paper the second approach was employed to 
derive a new set of effective potential force fields for 
molecules on a graphene sheet. The types of solid-fluid 
interactions of interest include the dispersion and repul- 
sive (Lennard- Jones, or LJ), the induction due to the 
dipole in the molecule and the polarizability of carbon 
atoms in graphene, and the interaction between multi- 
poles in the molecule and the permanent quadrupole in 
graphene. The potential functions derived in this work 
are extensions of a previous work by us j 1 61 ] for modeling 
interaction between gases and a semi infinite graphitic 
surface. The accuracy of the derived potentials are eval- 
uated by comparing the energies calculated from the inte- 
grated expressions and those from atomistic summation 
approach. The effectiveness of the potentials are demon- 
strated by using the derived formulas in simulating the 
adsorption of two representative gases into slits composed 
of single graphene sheets. 

In Section II, the functions are derived. Section III 
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and IV are examinations of the new potentials via com- 
parison with atomistic summation methods. Section V 
is concluding remarks. 



II. POTENTIAL DEVELOPMENT 

We assume that a molecule is located right above the 
center of a single graphene sheet. The origin of the sys- 
tem was set as the mass center of the molecule. If we 
only consider the pairwise interactions, the total poten- 
tial energy between the molecule and the graphene can 
be expressed by a summation 



U 



E 



u(ri 



(1) 



where = (xi,yi,z) denotes the position of a carbon 
atom in graphene relative to the molecule, uiti) is the 
pairwise potential between the molecule and carbon atom 
i in graphene, which is only dependent on the separation 
between them, rj = (xf + y 2 + z 2 ) 1 / 2 . For a graphene 
sheet of infinite size along its lateral directions (x and y), 
U only depends on the separation between the molecule 
and graphene along its surface normal z, and Eq.([T]) be- 
comes 



U[z) 



u(xi,yi,z). 



(2) 



For example, the Lennard- Jones interaction between the 
molecule and a carbon atom is given as 



u(r) = 4e sf 



(3) 



where e s f and tx s f are the LJ potential parameters for the 
cross interaction between the molecule and a graphene 
carbon atom. Substituting Eq.Q into Eq.([2]) gives 
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If we assume that the graphene sheet is homogeneous 
and continuous along its x and y directions, U(z) can be 
approximated by an integral over x and y 



U(z) = 4e sf dJ / - sf - * f dxdy, 

Jo Jo [{x 2 + y 2 + z 2 ) b (x^+y^+z 2 ) 6 

(5) 

where d s is the density of carbon atoms in a unit graphene 
area dxdy. The integration can be simplified by letting 
x 2 + y 2 = S, which gives dxdy = ndS. Then Eq.© 
reduces to 

U(z) = 47T£ sf d, 



(S + z 2 ) 6 (S + z 2 ) 3 



dS. (6) 



Eq. © leads to an integrated potential for the L J interac- 
tion between a molecule and an infinitely large graphene 
sheet 



U h j(z) = 2Tre si a 2 f d s 



(7) 



The derived expression Eq.© for the LJ interaction 
between a molecule and a graphene is analogous to a po- 
tential by Steele The interaction is only a function 
of the distance of the molecule from the graphene along 
its planar normal, z, thus the numerical calculation is 
substantially simplified. It is noteworthy to point out 
that unlike the potential for graphite [15[ the new func- 
tion for graphene, Eq.©, does not contain any empirical 
term that requires fitted parameters. This is due to the 
fact that no approximation was involved in the deriva- 
tion to account for the semi-infinite layers of graphene 
for a graphite surface. Only assumption here is that the 
atoms in graphene sheet are continuous, otherwise the 
derivation is strictly exact. 

For LJ particles, Eq. ([7]) should be a sufficient approx- 
imation for calculating the molecule-graphene potential. 
However, for molecules that has partial charges or per- 
manent multipole moments, additional polar interactions 
have to be included. For this purpose, similar procedures 
are applied to derive integrated forms for the induction 
and multipolar interactions between the molecule and 
a graphene sheet. We note that the above derivation 
procedure is applicable as long as the pairwise poten- 
tials are dependent only of r by an inverse power r~ n 
with n > 1. This condition is satisfied if we use point 
or angle-averaged dipoles and point quadrupoles for the 
molecule-graphene interaction. For example, the angle- 
averaged dipole-induced dipole interaction between a po- 
lar molecule and a carbon atom in the graphene is given 
by E3 



(4tc ) 2 



(8) 



where /if is the permanent dipole moment of the 
molecule, a s is the isotropic polarizability of a carbon 
atom in graphene, and eo is the vacuum permittivity. 

The angle- averaged dipole-quadrupole interaction is 
given by 



/cT(47re ) 2 r 8 



(9) 



where O s is the permanent quadrupole moment on each 
carbon atom in graphene, k is Boltzmann's constant and 
T is the absolute temperature. 

Finally, the angle-averaged quadrupole-quadrupole in- 
teraction for the molecule and carbon atom is given by 
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TABLE I: The Lennard- Jones and multipole parameters for 
the example molecules studied. 





£f (kcal/mol) 


at (A) 


/i f (Debye) 9 f (10" 20 C A 2 ) 


CH 4 


0.2378 


3.527 




H 2 


0.1554 


3.152 


2.351 


Cl 2 


0.7101 


4.115 


10.79 



where 0f is the permanent quadrupole moment of the 
molecule. 

Substituting Eqs. ([8"])-(jT0]l into Eq.© and integrat- 
ing, we obtain the following expressions. The integrated 
induction potential is 
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The integrated dipole-quadrupole potential is 
_ nd s tf 6 2 1 

The integrated quadrupole-quadrupole potential is 
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Likewise, the potential functions in Eq. (fTT| — (fl~3|) are 
only dependent on the separation between the multipole 
and graphene along its planar normal z. By using Eq.([7|) 
and Eqs.tfn Tj — (fl"3|) . one can calculate the interaction po- 
tential energy between a polar molecule and a graphene 
sheet without using the costly pairwise atomistic sum- 
mation. 



III. COMPARISON WITH ATOMISTIC 
POTENTIALS 

To examine the accuracy of the derived functions, we 
calculated the potential energies of several representa- 
tive nonpolar and polar molecules interacting with a sin- 
gle graphene sheet using the integrated functions Eqs.© 
and (fTTj) — (fl3|) . to compare with those from atomistic 
summation of potentials Eqs. ([3]) and (|8"j) — (fTU]) . The 
graphene sheet in the atomistic summation was modeled 
as a square about 10 nm in a side, containing 3680 car- 
bon atoms. Trial calculations indicate that such a system 
size contains more than enough carbon atoms to approx- 
imate a graphene sheet that is infinite along its lateral 
directions. A single molecule was placed over the center 
of the graphene sheet and its interaction potential en- 
ergies were computed by summing up interactions from 
each carbon atom. 

The LJ and multipole potential parameters for the ex- 
ample molecules chosen are shown in Table HI The nonpo- 
lar molecule studied is methane, modeled by a single LJ 




FIG. 1: Interaction of methane with a single graphene 
sheet calculated using the atom-explicit pairwise summation 
of Eq.([3]) (line) and the integrated function Eq.© (symbols). 



site [18| . The representative polar molecules selected are 
water and chlorine. Water is modeled by the SPC/E po- 
tential 19], with a strong dipole moment of 2.351 Debye. 
Please note that the LJ parameters of water in Tableware 
for the O atom. The chlorine molecule is modeled by a 
single LJ site plus a quadrupole moment of 10.79xl0~ 20 
C A 2 [20]. The purpose of the calculation was to study 
the accuracy of the integrated potential compared with 
atomistic summation. Therefore, the potential models 
for the molecules were selected arbitrarily. The values 
of potential parameters for carbon atoms in graphene 
are cx s =3.40 A, e s =0.05569 kcal/mol, 4=0.382 A" 2 [H], 
a c =1.76 A 3 [2l|, 6 s =-3.03xl0- 20 C A 2 ;22|. The cross 
interaction parameters £ s f and a s { were calculated using 
the Lorentz-Bertholet rules. 

In Figs. Q]to[3]we compare the potential energies calcu- 
lated from the integrated potentials and the atom-explicit 
potentials at 300 K, for CH4, H2O, and CI2. It is seen 
that the integrated expressions are in excellent agreement 
with the atomistic summation results. If we take the 
atomistic summation results as the standard, the typi- 
cal errors in the integrated potential are less than 2% at 
the energy minimum. The integrated functions always 
give a slightly more attractive potential compared to the 
atom-explicit methods, possibly due to the fact that we 
treat the graphene carbon atoms as continuum in the 
integrated potential. 

The integrated potentials are most accurate when the 
graphene sheet is infinitely large. However, we also car- 
ried out series of calculations using graphene sheets of 
various sizes, to test the applicability of the derived 
functions to the finite-size graphene. We found that 
the integrated potentials are in excellent agreement with 
the atomistic potentials as long as the square-shaped 
graphene sheet contains more than ~500 carbon atoms 
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FIG. 2: Interaction of a water molecule with a sin- 
gle graphene sheet calculated using the atomistic summa- 
tion (lines) and integrated functions (symbols). Atom- 
istic summations: LJ interaction Eq.@ (solid line), LJ 
plus induction interactions Eq. ([3J +Eq. ([8J (long dashed line), 
and LJ plus induction plus dipole-quadrupole interactions 
Eq.©+Eq.([SJ)+Eq.@ (dashed line). Integrated functions: 
LJ interaction Eq.© (circle), LJ plus induction interactions 
Eq.(Tr))+Ea. (|lll) (square), and LJ plus induction plus dipole- 
quadrupole interactions Ea.(l7)l+Ea. (fTT1l +Eq. (fT2)l (diamond). 




Z, A 

FIG. 3: Interaction a chlorine molecule with a single 
graphene sheet calculated using the atomistic summation 
(lines) and integrated functions (symbols). Atomistic summa- 
tions: LJ interaction Eq.© (solid line), LJ plus quadrupole- 
quadrupole interactions Eq . (f3|) +Eq ITOl (long dashed line). In- 
tegrated functions: LJ interaction Eq.([7|) (circle), LJ plus 
quadrupole-quadrupole interactions Eq . (fTt+Eq. (|13fl (square). 



(or a graphene sheet of ~38x38 A 2 ). This suggests that 
the derived potentials can also be readily used in esti- 
mating the interaction of molecules with finite graphene 
sheets that arc nanometer-scale in size. 

For the different types of interactions between a 
molecule and graphene, it is interesting to see that, for ex- 
ample from Fig. [5J the induction and dipole-quadrupole 
interactions for strongly polar molecules such as water 
are non-negligible compared with the LJ energy. On the 
other hand, based on the CI2 example, the quadrupole- 
quadrupole interaction is usually not significant, even for 
molecules with a strong quadrupole moment like chlorine. 



IV. APPLICATIONS IN ADSORPTION 
SIMULATION 

One of important applications of the potential func- 
tions developed in this work is molecular simulation 
of adsorption of molecules in pores composed of single 
graphene sheets. In turn, such simulations serve as a veri- 
fication of the accuracy and effectiveness of the integrated 
functions in practical applications. For this purpose, we 
chose to simulate the adsorption of two typical nonpolar 
and polar gases, methane and water, confined in single 
graphene sheets. The simulations were performed using 
the grand canonical Monte Carlo (GCMC) method, the 
detailed description of which can be found in literature 
23]. 

The potentials parameters describing the interactions 
of water and methane with graphene are given in Table HI 
The methane is modeled as a simple one-site LJ particle, 
with a cutoff of 9 A. The LJ interaction between water 
molecules is also modeled by a cutoff of 9 A, without 
long-range correction applied, as suggested by the origi- 
nal literature The electrostatic interaction between 
water molecules are modeled by the partial charges dis- 
tributed on the three charge sites on the water model, 
without long range correction. The interaction between 
the molecules and graphenes is calculated by either the 
atomistic summation approach or integrated functions 
for comparison. 

The GCMC cell for adsorption simulations is a rect- 
angular box. Two graphene sheets are placed in the cell 
with their planes parallel to the x — y plane of the box. 
The distance between the two graphene layers along the 
z direction is 15 A. Again, this separation is chosen ar- 
bitrarily. Periodic boundary conditions are applied in all 
three directions. Therefore, the height of simulation box 
in the z direction is 30 A. The typical box sizes in the 
x and y directions are set as 32 and 34 A respectively, 
leading to a cell volume of about 33 nm 3 . The types 
of move attempted during a GCMC simulation were se- 
lected randomly with probability of 0.40, 0.40, 0.10, and 
0.10 for displacements, rotations, creations, and deletions 
of a fluid molecule respectively. For simulation of the 
spherical methane, the rotation moves were merged to 
displacement moves. Each simulation included equilibra- 
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FIG. 4: Adsorption isotherm of methane in singlet graphene 
sheets at 77 K. Circles were calculated using the atomistic 
summation potential and squares were calculated using the 
integrated potential. Adsorbed amount (p*) and chemical 
potential (/x*) are in reduced unit. 

tion of 2xl0 6 MC moves and production of 2xl0 6 MC 
moves. 

The simulation of methane adsorption in graphene 
sheets were performed at 77 K [18| . Previous theoretical 
studies predict that water only wets graphitic carbon sur- 
faces when temperature is above ~510 K (UHII- There- 
fore, we chose to perform simulations of water/graphene 
at 550 K. In GCMC simulations the reduced chemical po- 
tential was varied to obtain isotherms and density profiles 
of adsorbed methane or water, using either atomistic or 
integrated potential for comparison. 

Calculated adsorption isotherms for methane in 
graphene sheets at 77 K are presented in Fig. [4] The 
circles are simulation results computed from the atom- 
istic summation of Eq.Q, while the squares are from 
simulations using the integrated potential Eq.©. It can 
be seen that, within statistical uncertainty, the isotherms 
from these two different potentials are in excellent agree- 
ment. The first plateau in the isotherm corresponds to 
the first layers of methane adsorbed on graphene sheets, 
and the second plateau is from the methane adsorbed 
onto first layers as chemical potential increases. This is 
confirmed by density profile distribution of the adsorbed 
phase shown in Fig. [5] The transition from nonadsorp- 
tion to the first layer occurs at about \x* — —8.6, and 
the transition from the first to the second layer occurs at 
about fi* = -7.3. 

It is known that the isotherm shapes, especially the 
layering transition of adsorbed phase calculated from 
simulations, are very sensitive to the solid-fluid interac- 
tions (26[. A small deviation in the solid- fluid potential 
can result in significant shift in the layer transition loca- 
tion in isotherms. The excellent agreement in the adsorp- 
tion layering transition for the two potentials indicates 



FIG. 5: Density profiles of methane adsorption in graphene 
sheets at 77 K and reduced chemical potential of —6.5. (a) 
was calculated using the atomistic summation potential and 
(b) was calculated using the integrated potential. Two curves 
are almost identical so (a) was shifted by 0.01 in P(z) for 
clarity. The dashed lines at z = ±7.5 A represent the two 
graphene sheets in the system. 
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FIG. 6: Adsorption isotherm of water in graphene sheets at 
550 K. Circles were calculated using the atomistic summation 
potential and squares were calculated using the integrated 
potential. Adsorbed amount (p*) and chemical potential (p*) 
are in reduced unit. 

that the integrated functions are robust alternatives to 
the atomistic models. 

One example of density profiles of methane adsorbed 
in graphene layers is given in Fig. [5] Sharp peaks are 
observed at z = ±4, ±11 A, which correspond to the 
first adsorbed methane layers on the graphene sheets. 
Additional broad peaks are observed at z = 0, ±15 A, 
representing the adsorbed second layers between the first 
layers. Again, it is seen that the density distribution 
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FIG. 7: Density profiles of water adsorption in graphene 
sheets at 550 K and reduced chemical potential of —26. (a) 
was calculated using the atomistic summation potential and 
(b) was calculated using the integrated potential. Two curves 
are almost identical so (a) was shifted by 0.005 in P(z) for 
clarity. The dashed lines at z = ±7.5 A represent the two 
graphene sheets in the system. 

of adsorbed methane in graphenes calculated using the 
atomistic and integrated potentials are in excellent agree- 
ment. Note that the density profiles calculated from the 
two potentials are almost identical to each other, so that 
we have to present them by shifting one curve along the 
vertical axis by a constant value to have a clear view. 

The results for water adsorption in singlet graphenes 
at 550 K are shown in Figs. [5] and [7] It is seen from 
the water adsorption isotherm that a continuous wet- 
ting transition occurs at about /x* = —33.4. The density 
profile peaks corresponding to the first adsorbed layers 
on the graphene sheets locate at z = ±4.4 and ±10.6 
A. Slightly different from that of methane, the adsorbed 
phase of water forms two second layers between the first 
layers, at z — ±1.5 and ±13.5 A, due to the relatively 
smaller size of water molecule compared to methane. The 
isotherms and density profiles of water in graphene sheets 
calculated from the integrated and atomistic potentials 
are also in excellent agreement. This indicates that the 
integrated potentials for multipolar interactions are also 
excellent approximations to the atomistic approach. 

We monitored the computational time required for the 
simulations using two different approaches. It is found 
that typically the computational time for the integrated 
potential is < 1% of that required for the atomistic sum- 
mations. 



V. CONCLUSIONS 

A set of effective potentials were derived for calcu- 
lating the Lennard-Jones, dipole-induced dipole, dipole- 
quadrupole, and quadrupole-quadrupole interactions be- 
tween fluid molecules and a graphene sheet. The in- 
tegrated potential functions depend only on the sepa- 
ration between a molecule and the graphene surface. 
They are mathematically simple and easy to use in ei- 
ther estimating the interaction between a single molecule 
and a graphene sheet or in large-scale molecular simula- 
tions. The potential energies calculated from the inte- 
grated potentials are in excellent agreement with the re- 
sults calculated from direct atomistic summations, while 
the integrated ones are computationally negligible com- 
pared with their atomistic counterparts. Adsorption sim- 
ulations of two representative gases in singlet graphene 
sheets were performed to further test the derived poten- 
tial functions. The adsorption isotherms and density pro- 
files calculated from the derived potentials and atomistic 
models are in excellent agreement, indicating that the 
derived potential can predict both the equilibrium and 
the structural properties of adsorbed phase with excellent 
accuracy. The potentials developed in this work provide 
a simple, accurate, and robust method for calculating 
the physical interaction of molecules and single graphene 
sheets. 

Finally, we note that the self-consistency in the polar- 
ization of carbon atoms by a polar fluid molecule was ne- 
glected in the derivation. However, we believe that the 
effect is relatively small compared with the dispersion, 
induction, and dipole-quadrupole interactions. Also, by 
using an integrated potential, one assumes the graphene 
surface is smooth and the impact of the surface corruga- 
tion on the molecule-solid interaction is ignored. 
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